function X = sample(mdl, n)
%SAMPLE Samples from the Gaussian model
%
% [ Syntax ]
%   - X = sample(mdl, n)
%
% [ Arguments ]
%   - mdl:          The Gaussian model
%   - n:            The number of vectors to sample
%   - X:            The matrix of sampled vectors
%
% [ Description ]
%   - X = sample(mdl, n) samples n vectors from the Gaussian model mdl.
%
% [ Remarks ]
%   - The implementation is based on randn and chol. If the covariance
%     matrix is semi-definite, it will use spectral decomposition.
%
% [ History ]
%   - Created by Dahua Lin, on Dec 23, 2007
%

%% parse and verify input

assert(isnumeric(n) && isscalar(n) && n >= 0, ...
    'n should be a non-negative integer scalar.');

%% main

if n > 0
    X = randn(getdim(mdl), n);
    
    [R, p] = chol(mdl.sigma);
    if p == 0   % full rank
        R = R';
    else        % rank deficient
        [evals, U] = slsymeig(mdl.sigma);
        evals(evals < 0) = 0;
        
        R = bsxfun(@times, U, sqrt(evals'));
    end
    
    X = R * X;
    
    if ~slisallzeros(mdl.mu)
        X = bsxfun(@plus, mdl.mu, X);
    end
    
else
    X = zeros(getdim(mdl), 0);
end